-- These routines use the "downhill simplex method". A simplex is a N dimensional geometrical figure with N+1 vertices. The algorithim evaluates the function at each vertex and steps the simplex toward the minimum by: 1) reflecting a vertex away from the function's highest point. 2) reflection and expansion away from the highest point. 3) contraction away from the highest point. 4) contraction in all dimensions toward the lowest point.
-- step until parameters change by less than tol
minimize = init,
step,
step while pchange > tol,
minp:=psum/m
pchange = sum((abs(p[lowest]-p[highest])/
(abs(p[lowest])+abs(p[highest])))[i],i,
1,n)
-- move simplex one step
step = rank,
try(-alpha), -- reflect from highest
try(gamma) when ytry ≤ y[lowest], -- keep going
(ysave:=y[higest],
try(beta), -- contract from highest
shrink when ytry ≥ ysave) when ytry ≥ y[high]
-- find lowest, highest and next highest vertex
rank = lowest:=1,high:=1,highest:=2,i:=1,
(lowest:=i when y[i]<y[lowest],
(high:=highest,highest:=i) when y[i]>y[highest],
high:=i when y[i]>y[high] and i ≠ highest,
i:=i+1) while i≤m
-- move highest vertex through face by fac
try(fac) = ptry:=psum*(1-fac)/n-
p[highest]*((1-fac)/n-fac),
ytry:=f(ptry),
(y:=replace(y,highest,ytry),
psum:=psum+ptry-p[highest],
p:=replace(p,highest,ptry)) when ytry<y[highest]
-- shrink entire simplex
shrinkp(k)[j] = .5*(p[k,j]+p[lowest,j]) dim[n]
shrink = kk:=1,
((psum:=shrinkp(kk),
p:=replace(p,kk,psum),
kk:=kk+1) when kk≠lowest) while kk≤m,
y:=fp
-- define finite arrays so they can be assigned
sump[j] = sum(p[i,j],i,1,m) dim[n]
initp[i,j]= guess[j]+(deltas[j] when i=j,0) dim[m,n]
fp[i]=f(p[i]) dim[m]
init = n:=count(guess), m:=n+1,
p:=initp,
psum:=sump,
y:=fp
replace(array,index,newval)[i] = newval when i=index,